# Set the working directory

setwd(wd)
source("~//DALY_function.R")
source("~/cost_functions.R")

#read CSV
list_file <- list.files()
list_markov <- list()
name_markov <- c()
for (i in c(1:length(list_file))){
  list_markov[[i]] <- read.csv(list_file[i], header = TRUE, sep = ",", quote = "",row.names("cycle"),
                                                              dec = ".", fill = TRUE, comment.char = "")
  name_markov[i] <- substr(list_file[i],1,(nchar(list_file[i])-8))
}
names(list_markov) <- name_markov


#get mean markov
list_markov_mean <- list_markov
for (i in c(1:length(list_markov))){
  list_markov_mean[[i]] <- get_mean_markov(list_markov[[i]])
}

#get Years of life saved
list_YLS <- list_markov_mean
for (i in c(1:length(list_YLS))){
  list_YLS[[i]] <- get_YLS(list_markov_mean$WHO_wVL, list_markov_mean[[i]])
}

#get Years of life saved, discounted
list_YLS_disc <- list_markov_mean
for (i in c(1:length(list_YLS_disc))){
  list_YLS_disc[[i]] <- get_YLS_disc(list_markov_mean$WHO_wVL, list_markov_mean[[i]])
}

# Create disability Weight index
col_names <- colnames(list_markov_mean[[1]])
DW_list_trt_gbd <- vector("list",length = length(col_names)-1)
names(DW_list_trt_gbd) <- col_names[1:(length(col_names)-1)]
for (i in c(1:length(DW_list_trt_gbd))) {
  DW_list_trt_gbd[i] <- 0
  if (grepl("\\.trt",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[i] <- 0
  if (grepl("CC.trt",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[i] <- 0
  if (grepl("DCC",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[i] <- 0.178
  if (grepl("HCC",names(DW_list_trt_gbd[i])) == TRUE) DW_list_trt_gbd[i] <- 0.288
}
DW_list_trt_shev <- DW_list_trt_gbd
for (i in c(1:length(DW_list_trt_shev))) {
  if (grepl("\\.ntrt",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0
  if (grepl("CC.trt",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0
  if (grepl("CC.ntrt",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0
  if (grepl("DCC.trt",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0.123
  if (grepl("DCC.ntrt",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0.250
  if (grepl("HCC",names(DW_list_trt_shev[i])) == TRUE) DW_list_trt_shev[i] <- 0.687
}
DW_list_trt_gbd_SA_1 <- DW_list_trt_gbd
for (i in c(1:length(DW_list_trt_gbd_SA_1))) {
  DW_list_trt_gbd_SA_1[i] <- 0
  if (grepl("\\.trt",names(DW_list_trt_gbd_SA_1[i])) == TRUE) DW_list_trt_gbd_SA_1[i] <- 0.025
  if (grepl("CC.trt",names(DW_list_trt_gbd_SA_1[i])) == TRUE) DW_list_trt_gbd_SA_1[i] <- 0.025
  if (grepl("DCC",names(DW_list_trt_gbd_SA_1[i])) == TRUE) DW_list_trt_gbd_SA_1[i] <- 0.178
  if (grepl("HCC",names(DW_list_trt_gbd_SA_1[i])) == TRUE) DW_list_trt_gbd_SA_1[i] <- 0.288
}

DW_list_trt_gbd_SA_2 <- DW_list_trt_gbd
for (i in c(1:length(DW_list_trt_gbd_SA_2))) {
  DW_list_trt_gbd_SA_2[i] <- 0
  if (grepl("\\.trt",names(DW_list_trt_gbd_SA_2[i])) == TRUE) DW_list_trt_gbd_SA_2[i] <- 0.05
  if (grepl("CC.trt",names(DW_list_trt_gbd_SA_2[i])) == TRUE) DW_list_trt_gbd_SA_2[i] <- 0.05
  if (grepl("DCC",names(DW_list_trt_gbd_SA_2[i])) == TRUE) DW_list_trt_gbd_SA_2[i] <- 0.178
  if (grepl("HCC",names(DW_list_trt_gbd_SA_2[i])) == TRUE) DW_list_trt_gbd_SA_2[i] <- 0.288
}

DW_list_trt_gbd_SA_3 <- DW_list_trt_gbd
for (i in c(1:length(DW_list_trt_gbd_SA_3))) {
  DW_list_trt_gbd_SA_3[i] <- 0
  if (grepl("\\.trt",names(DW_list_trt_gbd_SA_3[i])) == TRUE) DW_list_trt_gbd_SA_3[i] <- 0.075
  if (grepl("CC.trt",names(DW_list_trt_gbd_SA_3[i])) == TRUE) DW_list_trt_gbd_SA_3[i] <- 0.075
  if (grepl("DCC",names(DW_list_trt_gbd_SA_3[i])) == TRUE) DW_list_trt_gbd_SA_3[i] <- 0.178
  if (grepl("HCC",names(DW_list_trt_gbd_SA_3[i])) == TRUE) DW_list_trt_gbd_SA_3[i] <- 0.288
}

DW_list_trt_gbd_SA_4 <- DW_list_trt_gbd
for (i in c(1:length(DW_list_trt_gbd_SA_4))) {
  DW_list_trt_gbd_SA_4[i] <- 0
  if (grepl("\\.trt",names(DW_list_trt_gbd_SA_4[i])) == TRUE) DW_list_trt_gbd_SA_4[i] <- 0.1
  if (grepl("CC.trt",names(DW_list_trt_gbd_SA_4[i])) == TRUE) DW_list_trt_gbd_SA_4[i] <- 0.1
  if (grepl("DCC",names(DW_list_trt_gbd_SA_4[i])) == TRUE) DW_list_trt_gbd_SA_4[i] <- 0.178
  if (grepl("HCC",names(DW_list_trt_gbd_SA_4[i])) == TRUE) DW_list_trt_gbd_SA_4[i] <- 0.288
}

# Calculate cost

col_names <- colnames(list_markov_mean[[1]])
l_cost <- vector("list",length = length(col_names)-1)
names(l_cost) <- col_names[1:(length(col_names)-1)]
p_HCC_DCC_no_care <- 0.465
cost_DCC <- 251*p_HCC_DCC_no_care
cost_HCC <- 251*p_HCC_DCC_no_care

l_v_cost_NH <- list("out_ntrt_EP"=0,"out_ntrt_EN"= 0,"out_trt_EP"= 0, "out_trt_EN" = 0,
                    "Cirr"=167.4,"DCC"=251,"HCC"=251)
#TREAT-B : l_v_cost_no_VL
l_v_cost_no_VL <- list("out_ntrt_EP"=21.7,"out_ntrt_EN"= 14.2,"out_trt_EP"= 55.1, "out_trt_EN" = 55.1,
                       "Cirr"=167.4,"DCC"=251,"HCC"=251)
l_v_cost_WHO_VL <- list("out_ntrt_EP"=53,"out_ntrt_EN"= 45.5,"out_trt_EP"= 55.1, "out_trt_EN" = 55.1,
                        "Cirr"=167.4,"DCC"=251,"HCC"=251)
l_v_cost_EASL <- list("out_ntrt_EP"=59.1,"out_ntrt_EN"= 51.6,"out_trt_EP"= 55.1, "out_trt_EN" = 55.1,
                      "Cirr"=167.4,"DCC"=251,"HCC"=251)
l_v_cost_treat_all <- list("out_ntrt_EP"=0,"out_ntrt_EN"= 0,"out_trt_EP"= 55.1, "out_trt_EN" = 55.1,
                           "Cirr"=167.4,"DCC"=251,"HCC"=251)

l_cost_NH <- get_l_variable_cost(l_cost, l_v_cost_NH)
l_cost_no_VL <- get_l_variable_cost(l_cost,l_v_cost_no_VL)
l_cost_WHO_VL <- get_l_variable_cost(l_cost,l_v_cost_WHO_VL)
l_cost_EASL <- get_l_variable_cost(l_cost,l_v_cost_EASL)
l_cost_treat_all <- get_l_variable_cost(l_cost,l_v_cost_treat_all)


#cost detail
l_v_cost_detail_NH <- list("TDF"=0,"FU_EP_ntrt"= 0, "FU_EN_ntrt"= 0,"FU_EP_trt"= 0, "FU_EN_trt"= 0,"hospit_Cirr"=167.4,
                               "hospit_DCC"=cost_DCC,"HCC"=cost_HCC)
l_v_cost_detail_treatB <- list("TDF"=40.3,"FU_EP_ntrt"= 24.3, "FU_EN_ntrt"= 14.2,"FU_EP_trt"= 14.2, "FU_EN_trt"= 14.2, "hospit_Cirr"=167.4,
                               "hospit_DCC"=cost_DCC,"HCC"=cost_HCC)
l_v_cost_detail_WHO_VL <- list("TDF"=40.3,"FU_EP_ntrt"= 52.6, "FU_EN_ntrt"= 44.8,"FU_EP_trt"= 14.2, "FU_EN_trt"= 14.2, "hospit_Cirr"=167.4,
                               "hospit_DCC"=cost_DCC,"HCC"=cost_HCC)
l_v_cost_detail_EASL <- list("TDF"=40.3,"FU_EP_ntrt"= 53, "FU_EN_ntrt"= 51.5,"FU_EP_trt"= 14.2, "FU_EN_trt"= 14.2, "hospit_Cirr"=167.4,
                               "hospit_DCC"=cost_DCC,"HCC"=cost_HCC)
l_v_cost_detail_UTT <- list("TDF"=40.3,"FU_EP_ntrt"= 0, "FU_EN_ntrt"= 0,"FU_EP_trt"= 14.2, "FU_EN_trt"= 14.2, "hospit_Cirr"=167.4,
                             "hospit_DCC"=cost_DCC,"HCC"=cost_HCC)

cost_detail_NH <- get_cost_detail(l_v_cost_detail_NH,apply(list_markov_mean$`Natural History`,2,sum))
cost_detail_treatB <- get_cost_detail(l_v_cost_detail_treatB,apply(list_markov_mean$`TREATB_2`,2,sum))
cost_detail_WHO_VL <- get_cost_detail(l_v_cost_detail_WHO_VL,apply(list_markov_mean$`WHO_wVL`,2,sum))
cost_detail_EASL <- get_cost_detail(l_v_cost_detail_EASL,apply(list_markov_mean$`GS`,2,sum))
cost_detail_UTT <- get_cost_detail(l_v_cost_detail_UTT,apply(list_markov_mean$`Treat_all`,2,sum))

cost_detail <-  data.frame(cost_detail_NH,cost_detail_treatB,cost_detail_WHO_VL,cost_detail_EASL,
                           cost_detail_UTT)
#budget impact analysis using cost detail

cost_detail_NH_5y <- get_cost_detail(l_v_cost_detail_NH,apply(list_markov_mean$`Natural History`[c(1:5),],2,sum))
cost_detail_treatB_5y <- get_cost_detail(l_v_cost_detail_treatB,apply(list_markov_mean$`TREATB_2`[c(1:5),],2,sum))
cost_detail_WHO_VL_5y <- get_cost_detail(l_v_cost_detail_WHO_VL,apply(list_markov_mean$`WHO_wVL`[c(1:5),],2,sum))
cost_detail_EASL_5y <- get_cost_detail(l_v_cost_detail_EASL,apply(list_markov_mean$`GS`[c(1:5),],2,sum))
cost_detail_UTT_5y <- get_cost_detail(l_v_cost_detail_UTT,apply(list_markov_mean$`Treat_all`[c(1:5),],2,sum))

cost_detail_5y <-  data.frame(cost_detail_NH_5y,cost_detail_treatB_5y,cost_detail_WHO_VL_5y,cost_detail_EASL_5y,
                           cost_detail_UTT_5y)

#variable cost discount
NH_var_cost_disc <- get_var_cost_disc(list_markov_mean$`Natural History`,l_cost_NH)
GS_var_cost_disc <- get_var_cost_disc(list_markov_mean$GS,l_cost_EASL)
TreatB1_var_cost_disc <- get_var_cost_disc(list_markov_mean$TREATB_1,l_cost_no_VL)
TreatB2_var_cost_disc <- get_var_cost_disc(list_markov_mean$TREATB_2,l_cost_no_VL)
TreatB3_var_cost_disc <- get_var_cost_disc(list_markov_mean$TREATB_3,l_cost_no_VL)
TreatB4_var_cost_disc <- get_var_cost_disc(list_markov_mean$TREATB_4,l_cost_no_VL)
WHO_no_VL_var_cost_disc <- get_var_cost_disc(list_markov_mean$WHO_NoVL,l_cost_no_VL)
WHO_VL_var_cost_disc <- get_var_cost_disc(list_markov_mean$WHO_wVL,l_cost_WHO_VL)
Treat_all_var_cost_disc <- get_var_cost_disc(list_markov_mean$Treat_all,l_cost_treat_all)

t_var_cost_disc <-  rbind(NH_var_cost_disc,GS_var_cost_disc,TreatB1_var_cost_disc,TreatB2_var_cost_disc,
                          TreatB3_var_cost_disc,TreatB4_var_cost_disc, WHO_VL_var_cost_disc, Treat_all_var_cost_disc)


# Calculation of YLL

list_YLL <- list_markov_mean
for (i in c(1:length(list_YLL))){
  list_YLL[[i]] <- get_YLL(list_markov_mean[[i]],LE = 61,age_init = 20)
}

#calculate YLL, discounted
list_YLL_disc <- list_markov_mean
for (i in c(1:length(list_YLL_disc))){
  list_YLL_disc[[i]] <- get_YLL_disc(list_markov_mean[[i]],LE = 61,age_init = 20,0.03)
}

#Calculation of YLD
list_YLD_trt_gbd <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd))){
  list_YLD_trt_gbd[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_gbd)
}
list_YLD_trt_shev <- list_markov_mean
for (i in c(1:length(list_YLD_trt_shev))){
  list_YLD_trt_shev[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_shev)
}
list_YLD_trt_gbd_SA_1 <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_1))){
  list_YLD_trt_gbd_SA_1[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_1)
}
list_YLD_trt_gbd_SA_2 <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_2))){
  list_YLD_trt_gbd_SA_2[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_2)
}
list_YLD_trt_gbd_SA_3 <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_3))){
  list_YLD_trt_gbd_SA_3[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_3)
}
list_YLD_trt_gbd_SA_4 <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_4))){
  list_YLD_trt_gbd_SA_4[[i]] <- get_YLD(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_4)
}
#Calculation of YLD, discounted
list_YLD_trt_gbd_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_disc))){
  list_YLD_trt_gbd_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_gbd,disc = 0.03)
}
list_YLD_trt_shev_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_shev_disc))){
  list_YLD_trt_shev_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_shev,disc = 0.03)
}
list_YLD_trt_gbd_SA_1_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_1_disc))){
  list_YLD_trt_gbd_SA_1_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_1,disc = 0.03)
}
list_YLD_trt_gbd_SA_2_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_2_disc))){
  list_YLD_trt_gbd_SA_2_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_2,disc = 0.03)
}
list_YLD_trt_gbd_SA_3_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_3_disc))){
  list_YLD_trt_gbd_SA_3_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_3,disc = 0.03)
}
list_YLD_trt_gbd_SA_4_disc <- list_markov_mean
for (i in c(1:length(list_YLD_trt_gbd_SA_4_disc))){
  list_YLD_trt_gbd_SA_4_disc[[i]] <- get_YLD_disc(list_markov_mean[[i]],DW_list = DW_list_trt_gbd_SA_4,disc = 0.03)
}
# calculation of DALY
list_DALY_trt_shev <- list_markov_mean
for (i in c(1:length(list_DALY_trt_shev))){
  list_DALY_trt_shev[[i]] <- list_YLD_trt_shev[[i]] + list_YLL[[i]]
}

list_DALY_trt_gbd <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd))){
  list_DALY_trt_gbd[[i]] <- list_YLD_trt_gbd[[i]] + list_YLL[[i]]
}

list_DALY_trt_gbd_SA_1 <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_1))){
  list_DALY_trt_gbd_SA_1[[i]] <- list_YLD_trt_gbd_SA_1[[i]] + list_YLL[[i]]
}

list_DALY_trt_gbd_SA_2 <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_2))){
  list_DALY_trt_gbd_SA_2[[i]] <- list_YLD_trt_gbd_SA_2[[i]] + list_YLL[[i]]
}

list_DALY_trt_gbd_SA_3 <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_3))){
  list_DALY_trt_gbd_SA_3[[i]] <- list_YLD_trt_gbd_SA_3[[i]] + list_YLL[[i]]
}

list_DALY_trt_gbd_SA_4 <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_4))){
  list_DALY_trt_gbd_SA_4[[i]] <- list_YLD_trt_gbd_SA_4[[i]] + list_YLL[[i]]
}

# calculation of DALYs, discounted
list_DALY_trt_shev_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_shev_disc))){
  list_DALY_trt_shev_disc[[i]] <- list_YLD_trt_shev_disc[[i]] + list_YLL_disc[[i]]
}

list_DALY_trt_gbd_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_disc))){
  list_DALY_trt_gbd_disc[[i]] <- list_YLD_trt_gbd_disc[[i]] + list_YLL_disc[[i]]
}

list_DALY_trt_gbd_SA_1_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_1_disc))){
  list_DALY_trt_gbd_SA_1_disc[[i]] <- list_YLD_trt_gbd_SA_1_disc[[i]] + list_YLL_disc[[i]]
}
list_DALY_trt_gbd_SA_2_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_2_disc))){
  list_DALY_trt_gbd_SA_2_disc[[i]] <- list_YLD_trt_gbd_SA_2_disc[[i]] + list_YLL_disc[[i]]
}
list_DALY_trt_gbd_SA_3_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_3_disc))){
  list_DALY_trt_gbd_SA_3_disc[[i]] <- list_YLD_trt_gbd_SA_3_disc[[i]] + list_YLL_disc[[i]]
}
list_DALY_trt_gbd_SA_4_disc <- list_markov_mean
for (i in c(1:length(list_DALY_trt_gbd_SA_4_disc))){
  list_DALY_trt_gbd_SA_4_disc[[i]] <- list_YLD_trt_gbd_SA_4_disc[[i]] + list_YLL_disc[[i]]
}

